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The development and verification of a one-dimensional material thermal response code 
with ablation is presented. The implicit time integrator, control volume finite element 
spatial discretization, and Newton’s method for nonlinear iteration on the entire system 
of residual equations have been implemented and verified for the thermochemical ablation 
of internally decomposing materials. This study is a continuation of the work presented 
in "One-Dimensional Ablation with Pyrolysis Gas Flow Using a Full Newton’s Method 
and Finite Control Volume Procedure" (AIAA-2006-2910), which described the derivation, 
implementation, and verification of the constant density solid energy equation terms and 
boundary conditions. The present study extends the model to decomposing materials 
including decomposition kinetics, pyrolysis gas flow through the porous char layer, and 
a mixture (solid and gas) energy equation. Verification results are presented for the 
thermochemical ablation of a carbon-phenolic ablator which invloves the solution of the 
entire system of governing equations. 


N omenclat ur e 


A Area 

A Area vector 

a , b Constants 

B' c Dimensionless char ablation rate 
Ch Corrected heat transfer Stanton number 
C ho Uncorrected heat transfer Stanton number 
C m Mass transfer Stanton number 
C p Specific heat at constant pressure 
C v Specific heat at constant volume 
DE Discretization error 
E Element energy convection rate 

E Energy content 
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e Specific internal energy 
(e) Element 
h Specific enthalpy 
k Thermal conductivity 
L Time dependent domain length 
to" Mass loss rate per unit area 
N Shape function 
Pr Prandtl number 
Q* Heat of ablation 
Q Element heat conduction rate 
q" Heat flux 
q" Heat flux vector 
R Residual 
r Recovery factor 
s Surface recession 
s Surface recession rate 
T Temperature 
t Time 

u Control volume boundary velocity 
u Velocity vector 
V Volume 

x Coordinate with respect to instantaneous ablation front 
z Coordinate with respect to initial ablation front 
z Nodal velocity 
a Thermal diffusivity 
S Perturbation 
e Emissivity 

Landau coordinate 
9 Transformed temperature 
A Blowing reduction parameter 
£ Local coordinate 
p Density 

a Stefan-Boltzman constant 
Stanton number correction 

Subscript 

abl Ablation 

ah Aerodynamic heating 

back Back boundary 

blw Blowing 

c Char 

cw Cold wall 

e Boundary layer edge condition 
full Full sensitivity matrix 
hw Hot wall 

iter Iterative solution procedure 
j Nodal index 

lag Property lagging solution procedure 

o Initial value 

r Recovery value 

rad Radiation 

ref Reference value 

res Reservoir property 

surf Surface 

tri Tridiagonal sensitivity matrix 
w Property of gases adjacent to the wall 
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Superscript 
n Time level 

v Iteration level 

* Reference value 


I. Introduction 

F OR nearly half of a century, numerical heat transfer and ablation modeling have been important tools 
in the design and analysis of rocket nozzles and re-entry vehicle heat-shields since they provide a means 
of simulating material thermal response including transient heat conduction, shape change, and in-depth 
decomposition. Accurate prediction of mass loss and energy transfer helps minimize vehicle weight while 
still maintaining the nozzle lining’s or heat-shield’s protective capabilities. In the 1960s, Moyer and Rinclal 1 
developed a one-dimensional, finite difference ablation code (Charring Materials Ablation, CMA) that em- 
ployed a translating grid scheme in which the grid was attached to the receding surface and the overall 
number of nodes in the domain would be reduced as mass was removed at the ablating surface. With this 
method, each node translates with the same velocity as the surface and the size of grid cells remains fixed 
except for the back boundary cell which reduces in size until the cell has reached a critical thickness. This 
"node-dropping" translating grid scheme was a common approach during the 1960s and 1970s and is still 
used today. 

Moyer and Rindal used a fully implicit time integrator and lagged the temperature dependent material 
properties and surface recession rate at the interior nodes while the nonlinear surface energy balance was 
solved iteratively to give the updated surface recession rate. Explicit integration of the solid energy equation 
has also been attempted by Moyer and Rindal during the development of CMA and more recently in a 
finite volume code by Suzuki et al. 2 In 1988, Blackwell 5 presented the control volume finite element 
method (CVFEM) for solving ablation problems on translating/node-dropping grids with a fully implicit 
time integrator with which he implemented the Moyer and Rindal approach to iterating only on the nonlinear 
surface energy balance. 

Later, a contracting grid scheme for one-dimensional ablation was presented by Blackwell and Hogan 3 4 
in which the relative size of grid cells and the total number of cells remains fixed throughout the problem, 
and each node translates at a fraction of the surface recession rate. This eliminates the node-dropping 
complexity, but adds a change in energy storage due to cell volume reduction. This method also has a more 
logical extension to multidimensional space as presented by Hogan et al. 5 Blackwell’s and Hogan’s work 
also implemented the CVFEM with a fully implicit time integrator that again iterated only on the surface 
energy balance to determine recession rate. More recent studies on one-dinrensional and multidimensional 
ablation modeling and application work can be attributed to Suzuki et al. 6 and Chen and Milos.' ,8,9 

The current study is a presentation of the development and verification of a one-dimensional thermal 
response and ablation code using the CVFEM, a contracting grid scheme, and fully implicit time integration. 
Although the model has been developed for planar, cylindrical, and spherical geometries, only the solution for 
planar geometries is presented. This is a simplification to show the development, verification, and usability of 
the method. The governing equations are solved using the block Gauss-Seidel segregated solution procedure. 
Newton’s method is used to iterate on the entire system of nodal residual equations to determine the time 
accurate temperature and density dependent properties and surface recession rate. Newton’s method is only 
used to solve the segregated mixture energy and gas phase continuity equations while direct integration of 
the decomposition kinetics is used to update the solid density. 

II. Governing Equations 

The equations that govern the solid/gas system of the porous charring ablator include solid energy and 
continuity equations as well as the Navier-Stokes equations as applied to all of the gaseous species considered. 
In the general case, it is possible that the pyrolysis gases react among themselves, erode the remaining solid, 
or deposit residue (coke) on the solid, but these phenomena are neglected. Under the assumptions that 
the pyrolysis gas is a single nonreactive entity, the solid and gas are in thermal equilibrium, and there is 
no in-depth energy source, then the solid and gas energy equations for a moving grid reduce to a mixture 
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energy equation given by 


q" • dA + I cj) p g h g v g • dA — I phv cs • dA + — / pedV = 0 


dt 


conduction gas flux 

and the mixture continuity equation is 


grid convection 


energy content 


(t)p g Vg - dA - I pv cs • dA+ 


dt 


pdV =0 


gas flux 


grid convection mass content 


which is the sum of the solid and gas phase continuity equations given by 



solid mass content grid convection solid mass source 


(i) 


(2) 


( 3 ) 


and 


gas: 



gas mass content 



gas flux grid convection gas mass source 


( 4 ) 


where 


m”'dV + 


( 5 ) 


Applying the Navier-Stokes momentum equations to flow through the char layer would require detailed 
knowledge of the pore structure, and that information is not known. Instead we use Darcy’s law which 
relates the volumetric flow rate (Q) of a laminar flowing fluid to the local pressure gradient within a fully 
saturated porous medium by 


Q = -A— VP 
d 


( 6 ) 


The superficial or filtration velocity is the volumetric flow rate averaged over the cross-sectional area of the 
medium and is given by 


v 


/ 


Q 

A 


— -VP 


( 7 ) 


Since the momentum equation for the pyrolysis gases is assumed to be Darcy’s law, which provides a 
simple relationship between fluid velocity and the pressure gradient, Darcy’s law in conjunction with the 
perfect gas law is treated as a closure relationship for the energy-continuity system of equations with surface 
recession rate, temperature, solid density, and gas density being the dependent variables. As a result, there 
is no independent solution of a momentum equation since this is inherent in the gas phase continuity equation 
solution. 

The block Gauss-Seidel procedure uses Newton’s method for a system of nonlinear equations to solve 
the segregated mixture energy and gas phase continuity equations. Consequently, a residual formulation, 
or A-formulation, of the control volume energy and gas mass balance equations is solved resulting in the 
linear system’s Jacobian matrix, [A] in Eq. 8, being comprised of the residuals’ sensitivities to the dependent 
variables. 

[A] Ax = R (8) 


As a result, convergence can be monitored as the correction vector, Ax, or the residual vector, R, approaches 
zero. The solid continuity equation is not solved in this manner and is instead solved through direct 
integration of the decomposition kinetics. 

The solution procedure for the entire system of governing equations can be outlined as follows: 


1. Iteratively solve the mixture energy equation with nodal temperatures and surface recession rate as 
the dependent variables 
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• Adjust grid throughout iteration procedure (the contracting grid scheme is described in the pre- 
vious paper 10 and Amar 11 ) 

• Hold nodal density and gas velocity values constant 

2. Integrate kinetic relationships to solve solid phase continuity equation 

• Use updated grid and temperature held 

3. Iteratively solve the gas phase continuity equation with nodal gas densities as the dependent variables 

• Use Darcy’s law to determine pyrolysis gas velocity in the pore space 

• Hold solid density, temperature, and recession rate constant from previous converged solution of 
the energy and solid continuity equations 

This segregated solution procedure in which updated information is being used as soon as it is available 
is commonly referred to as the block Gauss-Seidel method which converges linearly if implemented correctly. 
Steps 1-3 are repeated until global convergence of the Gauss-Seidel method is reached, and the system is 
considered globally converged when no local nonlinear iteration is necessary to converge the mixture energy 
equation. 


III. Mixture Energy Equation 

derivation details will be in full paper 

IV. Solid Phase Continuity Equation 

Since the solid density gradient in the decomposition region is known to be much larger than the tem- 
perature gradients in the domain, a finer mesh is necessary for the solution of the solid phase continuity 
equation than was needed to solve the energy equation. CMA and several other studies employ a sub-mesh 
scheme in which there are a given number of nodelets (or sub-nodes) within each element on which the solid 
phase continuity equation is solved. Integral effects of the nodelet solutions are applied at nodes during the 
solution of the mixture energy equation. This method saves a significant amount of computational time 
since it avoids excessive grid points during the mixture energy equation solution. The nodelet scheme has 
not been implemented in the current model, therefore each governing equation is solved on the same mesh. 

Implicit integration of the solid mass conservation equation through the CVFEM and nodal mass balance 
equations is not performed. Instead, direct integration of the kinetic relationships accounting for the moving 
grid, which has been used in previous studies, 1 was chosen to determine the solution. The development 
of the fixed grid solid continuity equation solution procedure is presented, and the moving grid effects are 
considered afterwards. 

A. Solution Procedure 

The solid phase continuity equation on a fixed grid is given by 

f m'gdV 

J CV ^ 

solid mass source 

where 

Pa — T (. Pa + Pb) + (1 — r) p c 

' v ' ' v ' 

resin binder 

Substituting Eq. 10 into Eq. 9 and knowing that the component source terms can be volumetrically weighted 
like the component densities gives 

^ f [r (p A + Pb) + (1 - r) Pc \ dV= [ [T « + m"') + (1 - T) m£'] dV (11) 

^ J cv J CV 


( 9 ) 

(10) 



solid mass content 
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Since the decomposition of each component is independent of the other components when coupling effects 
from the other governing equations are ignored, the solid phase continuity equation can be further reduced 
to component continuity equations given by 


d 

dt 



rh'-'dV for i = A, B, and C 


(12) 


Now there are four equations (Eqs. 11 and 12) with four unknowns ( p A , p B , p c , and p s ), and the procedure 
is to solve the component continuity equations and then determine the solid density as a post-processing 
step via Eq. 10. The analogous finite difference relationship to Eq. 12 is 


dpj 

dt 



e ~Ei/RT 


for i = A, B, and C 


(13) 


Since the density history at a given location only depends on the temperature history at that location, Eq. 
13 can be integrated on a location-by-location basis (accounting for the moving grid) and the updated solid 
density can be determined. 


B. Integration of Kinetic Equations 

be rearranged to give 
B, and C 

where the dimensionless relative density is given by 

Pi Pet 

Wi = 

P Vi 


Since p v and p c . are constants, the kinetic equation in Eq. 13 can 

^ = -kiwf i e~ Ei / RT for * = A, 


(14) 


(15) 


Direct integration of Eq. 14 over one time step while treating temperature explicitly in time, gives 
Ui (^+\t" +1 ) = {[wi (1 

for 


Wi 


(16) 


and 


(*? +1 ,i n+1 ) = [wi (z™+\t n )]e 




-E i /R'r( Ij "+ 1 ,t") 


(17) 


for i/) i = 1 

Eqs. 16 and 17 were implemented in CMA to integrate the kinetic equations regardless of the temperature 
change over the time interval. They are applied in the current model only if the temperature remains 
constant over the time interval [T (z" +1 , t n+1 ) = T (zj +1 ,t n )] . 

An alternative approach to determining the dimensionless relative densities at time level n + 1 is to 
assume that the temperature rise rate over the time interval A t n+1 is a constant, 0. 


dT 

dt 


= 0 = 


2 


T(zJ +1 ,t n+1 ) -T(zJ +1 ,t n ) 
At 


(18) 


This is consistent with the implicit time integrator’s ability to at best capture a linear temperature rise rate 
since the method is first order. Using Eq. 14 and the chain rule gives 


dwi 

~dT 


0 = —kiwf^e E d RT f 01 - i = A, B, and C 


(19) 
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Integrating and solving for Wi (zj +1 ,t n+1 ) gives 






V’j — 1 


0 


/r(*ynt») 


e~ E ^ RT dT 


(20) 


and 


for 


( b. r T(zJ+\t”+ l ) ^ 

iH (4 +1 ,t n+1 ) = exp In [ Wi (+++") 1 - £ / e~ E */ RT dT 

V eJT(z 7 +\t») J 


(21) 


for ipi = 1 

/' T ( z ? +1 > t " +1 ) 

/ e~ Ei ^ RT dT can be determined through the use of exponential integrals or numerical integration 

techniques, and Simpson’s rule was chosen for this study. Since the new dimensionless relative densities 
are known, the component densities can be determined via Eq. 15, and the solid bulk density can now be 
determined from 

p n s +1 = r (pa +1 + Pb +1 ) + (i - E) Pc +1 (22) 

or in terms of dimensionless relative densities 


p : +1 = p cs +r(p. 


v a W A 


n+1 


+ Pvb W B +1 ) + (! - r ) Pvc W C 


n+1 


(23) 


Since problems are solved on contracting grids for which the spatial coordinate of each node changes while 
the Landau coordinate of each node is fixed, it is necessary to determine temperature and dimensionless rel- 
ative density profiles at the previous time level and current nodal locations (T (+ + , f") and vjj (+++")). 
Since properties are assumed to vary linearly within elements, linear interpolation in the old (time level n) 
profiles is used to calculate these properties. 

V. Gas Phase Continuity Equation 

derivation details will be in full paper 


VI. Verification Results: Thermo chemical Ablation of a Decomposing 

Material 


A. Problem Statement 

Since most of the energy equation terms and boundary conditions have been previously verified (see the 
previous paper 10 and Amar 11 ), the decomposing material thermochemical ablation problem is intended 
to verify the implementation of the gas phase continuity equation in addition to the pyrolysis gas effects 
in the mixture energy equation. Due to the difficulty in determining an analytical solution, Richardson 
extrapolation is used to approximate the exact solution, and the global convergence rate and nonlinear 
convergence rates of the mixture energy and gas phase continuity equations are determined. 

Consider a 0.5 inch thick planar slab of carbon-phenolic (with material properties given in Amar 11 ) that 
is subject to a typical heating load characteristic of a ballistic splrere-cone re-entry body. The Stanton 
number is corrected for both hot wall and blowing effects and the surface is subject to a far-field radiation 
condition where the source temperature is 414.0 °R. The back boundary is adiabatic and impermeable, and 
the time dependent specified surface pressure is provided by the aerodynamic heating boundary condition 
input. 

B. Grid Refinement Studies 

The domain was discretized with a series of four grids with At/Az% = constant, and the grid parameters 
can be seen in Table 1. Since the solid continuity equation is solved through direct integration of the 
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Table 1. Parameters for carbon-phenolic thermochemical ablation problem grid refinement study. 


Grid 

# elem 

A z Q , in 

At, s 

coarse 

50 

0.01 

0.1 

medium 

100 

0.005 

0.025 

fine 

200 

0.0025 

0.00625 

extra-fine 

400 

0.00125 

0.0015625 


decomposition kinetics, it is unclear how its solution and consequently its coupling with the other governing 
equations will behave as the mesh is refined. As a result, an alternative procedure to the previous verification 
problems was used in order to isolate the mixture energy and gas phase continuity equations from the solid 
phase continuity equation solution. The procedure is as follows: 

1. Solve the problem on the extra- fine mesh until a problem time of t = 25 s is reached 

2. The conditions found at t = 25 s can then be used to re-initialize the problem on the entire series of 
grids resulting in 

• Nonuniform initial profiles for: 

a. Temperature 

b. Solid density 

c. Porosity 

d. Permeability 

e. Gas density 

f. Gas generation rate 

g. Gas velocity 

• Contracted grid 

• Re-initialization during an ablating interval 

3. Verify the mixture energy equation 

• Integrate over 0.1 s (restarting at t = 25 s) for the entire series of grids only advancing the mixture 
energy equation in time while holding all other terms constant 

a. Integration over a small time interval reduces errors associated with inconsistent nodal locations 

b. Result is as revealing as integrating over a long time interval if all of the relevant terms are 

being employed which is ensured by the nonuniform initialization 

• Perform Richardson extrapolation at t = 25.1s on the surface temperature, recession rate, and 
back face temperature. The results showing the second order grid convergence can be seen in 
Figure 1. 
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Figure 1. Mixture energy equation grid convergence results for carbon-phenolic thermochemical ablation 
verification problem. 


4. Verify the gas phase continuity equation 

• Integrate over 1.0 s (restarting at t = 25 s) for the entire series of grids only advancing the gas 
phase continuity equation in time while holding all other terms constant 

• Perform Richardson extrapolation at t = 26.0 s on the surface gas flux, an interior gas density 
value at a consistent spatial location, and back face gas density. The results showing the second 
order grid convergence can be seen in Figure 2. 

a. An interior location is chosen for Richardson extrapolation on the gas density because the 

surface gas density is specified through the boundary conditions, z = 1.07 x 10~ 2 in is the 
position of the first interior node on the coarsest mesh which has a collocated node on each 
of the finer meshes. 

b. Integration over a longer time interval was performed since the mesh does not contract during 

the solution of the gas phase continuity equation 


C. Nonlinear Convergence Studies 

For the mixture energy equation, the error metric used to determine the nonlinear convergence rate is the 
recession rate as given by 

(error)" = (s)" - (, s) converged (24) 

and the error metric for the gas phase continuity equation is the relative L 2 norm of the density correction 
vector. 

( \v / r \ v it \ converged ir>r\ 

(error) = (L 2 ) - (L 2 ) (25) 

The convergence rates are shown in Figures 3 and 4, and the results follow the second order line. 

The block Gauss-Seidel method used to solve the entire system of governing equations is expected to 
exhibit linear convergence. Since the solution procedure achieves coupling through iteration, information 
can only be used as soon as it is updated and a smooth convergence curve should not be expected. To 
examine the global convergence, it is necessary to view the problem in global iteration space as opposed to 
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Number of Elements 


Figure 2. Gas continuity equation grid convergence results for carbon-phenolic thermochemical ablation ver- 
ification problem. 



Figure 3. Nonlinear convergence of surface recession rate for carbon-phenolic thermochemical ablation verifi- 
cation problem. 
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Figure 4. Nonlinear convergence of gas continuity equation’s L 2 error norm for the carbon-phenolic thermo- 
chemical ablation verification problem. 
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Figure 5. Global convergence rate for surface recession rate and gas density where v is the global iteration 
index. 


local (Newton) iteration space. Let the iteration index (u) now represent the global iteration. Since the 
solution globally converges to the exact solution plus some discretization error, the difference between the 
locally converged solution at a given global iteration (X v ) and the globally converged solution ( X 9 ' c ■) is used 
as the error metric. 

(error) 1 ' = X" - X g - C - (26) 

The global convergence rates of the surface recession rate and gas density are shown in Figure 5 while the 
global convergence rates for the surface temperature and solid density are shown in Figure 6. It is important 
to note that the surface gas density only changes with surface temperature since the boundary pressure is 
known from the boundary conditions. As a result, the surface gas density and temperature will always 
converge to machine precision in the same number of global iterations. On the other hand, the surface solid 
density reaches global convergence in fewer iterations and experience has shown that this is generally the 
case. 

D. Code-to-Code Comparison 

Although code-to-code comparisons are not part of the formal verification process, comparisons with estab- 
lished codes were performed to see how the improved scheme changes the final solution. The carbon-phenolic 
thermochemical ablation problem was solved on the medium grid in Table 1 using both CMA and the code 
developed during this study. Although CMA has the capability to integrate the decomposition kinetics on 
a finer mesh than is used to solve the mixture energy equation, this option was not exercised in order to 
make a more direct comparison. It is also important to note that the primary differences between the two 
codes for this simulation are as follows: 

1. CMA iterates only on the surface energy balance and lags temperature dependent properties and surface 
recession rate for the interior node solution while the research code iterates on the entire system of 
residual equations to determine time accurate nonlinearities 

2. CMA uses a translating/node-dropping grid scheme, and the research code uses a contracting grid 
scheme 

3. CMA uses a lumped capacitance method, and the research code uses a distributed capacitance method 
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Figure 6. Global convergence rate for surface temperature and solid density where v is the global iteration 
index. 
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Figure 7. Comparison of surface temperature predictions for the carbon-phenolic thermochemical ablation 
verification problem. 


4. CMA linearly interpolates on In ( B ' c ), and the research code linearly interpolates on B' c 

5. CMA assumes that all gas generated further in-depth than a given node passes that node in the time 
interval in which the gas was generated and the flow work is accounted for, and the research code 
allows internal pressure gradients to drive the flow of pyrolysis gases through the pore space 

6. CMA does not account for the gas phase internal energy in the energy capacitance term, and the 
research code does account for the gas phase internal energy 

The comparisons for the predicted surface temperature, back face temperature, recession, and recession 
rate histories can be seen in Figures 7, 8, 9, and 10 respectively. Due to the number of numerical differences 
between the codes and the number of physical phenomenon that are being modeled throughout the solution 
procedure, it is difficult to discern what modeling differences between the two codes is causing the discrepancy 
in the solutions. From these results it is evident that CMA’s pyrolysis gas assumptions, described in items 
5 and 6 in the above list, are valid since including the porous flow solution (as is done in the research code) 
does not significantly change the mixture energy equation solution. This suggests a dominantly one-way 
coupling between the mixture energy equation and the gas phase continuity equation, which is expected 
since gas velocities are low and the heat capacity and thermal conductivity for gases are generally much 
lower than they are for a solid. 

E. Porous Flow Results 

One of the primary goals of this study is the implementation and verification of the gas phase continuity 
equation with porous flow assumptions in order to predict in-depth pressure in a charring ablator. Figures 11- 
19 show profile histories for several properties in the charring carbon-phenolic ablator including gas density, 
pressure, and gas superficial velocity, which are predicted through solving the gas phase continuity equation 
with porous flow assumptions. It is evident that the gas density, pressure, and solid decomposition (or 
gas generation) rate profiles have steep gradients when compared to the temperature profile. This suggests 
that if a sub-mesh scheme is implemented, then both the solid and gas phase continuity equations should be 
solved on the sub-mesh. In addition, the apparent "kinks" in the permeability profiles in Figure 19 are a 
result of the available permeability data that was used for the simulation. 
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Time (s) 

Figure 8. Comparison of back face temperature predictions for the carbon-phenolic thermochemical ablation 
verification problem. 



Figure 9. Comparison of surface recession predictions for the carbon-phenolic thermochemical ablation veri- 
fication problem. 
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Figure 10. 
verification 


Figure 11. 



Comparison of surface recession rate predictions for the carbon-phenolic thermochemical ablation 
problem. 



Temperature profile history for the carbon-phenolic thermochemical ablation verification problem. 
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Figure 12. 


Figure 13. 



Solid density profile history for the carbon-phenolic thermochemical ablation verification problem. 



z (in) 


Gas density profile history for the carbon-phenolic thermochemical ablation verification problem. 


17 of 21 


American Institute of Aeronautics and Astronautics Paper 2007-???? 



Z (in) 


Figure 14. Pressure profile history for the carbon-phenolic thermochemical ablation verification problem. 
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Figure 15. Extent of reaction (/ 3 ) profile history for the carbon-phenolic thermochemical ablation verification 
problem. 
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Figure 16. Solid decomposition rate profile history for the carbon-phenolic thermochemical ablation verification 
problem. 



Figure 17. Gas superficial velocity profile history for the carbon-phenolic thermochemical ablation problem 
where a negative velocity is toward the ablating surface. 
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Figure 18. Porosity profile history for the carbon-phenolic thermochemical ablation verification problem. 
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Figure 19. Permeability profile history for the carbon-phenolic thermochemical ablation verification problem. 
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Figure 15 shows that ablation begins to occur before the surface solid becomes fully charred meaning that 
there is still non-decomposed resin in the solid. This creates a discrepancy since the surface thermochemical 
solution assumes that the solid participating in the reactions is fully charred. According to the decomposition 
kinetics, it takes an infinite amount of time for the composite to decompose to its char state. As a result, to 
ensure that the thermochemistry data is accurately describing the surface ablation phenomenon, a tolerance 
in the vicinity of the char density could be implemented with which the solid density is automatically adjusted 
to the char value once this tolerance has been reached. This method is implemented in CMA, but it adds 
the additional complication of accounting for this mass when solving the gas phase continuity equation. 

VII. Conclusions 

A one-dimensional charring ablator thermal response code with a contracting grid scheme has been 
developed and verified to exhibit second order spatial accuracy for the thermochemical ablation problem. 
Newton’s method for the entire system of equations has been implemented and has also been verified to 
exhibit second order nonlinear convergence rates for both the mixture energy and gas phase continuity 
equations. The model has shown good agreement with CMA, but it also includes several improvements in 
the solution procedure, which include the calculation of time accurate nonlinearities, improved integration 
of decomposition kinetics, and the prediction of pyrolysis gas flow and pore pressure. 
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